traits
above- and belowground
acquisitive vs. conservative resource-use strategies
collaborative vs. individual resource-use strategies
fine root proportion
thickness of roots
Goal: Examine how traits shift in community compositional (ordinational) space, placing each trait on an axis.
Author: Caryn D. Iwanaga
Updated: 01/24/2025
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.2 ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-8
cover.dat.labels <- read.csv("https://www.dropbox.com/scl/fi/up8nnzkcpchsr45f8cm92/Compost_Cover_LongClean.csv?rlkey=z2tvaj8t6khadef7ydz782zka&st=qwef9ys0&dl=1") %>%
mutate(
nut_trt = factor(ifelse(nut_trt =="C", "+compost", ifelse(nut_trt =="F", "+N fertilizer", ifelse(nut_trt =="N", "control", NA))), levels = c("control", "+N fertilizer", "+compost")),
ppt_trt = ifelse(ppt_trt == "D", "dry", ifelse(ppt_trt == "XC", "control", ifelse(ppt_trt == "W", "wet", NA))))
# Define a reusable color scale
# custom_colors <- scale_fill_manual(
# values = c(
# "control" = rgb(195, 197, 193, maxColorValue = 250),
# "+N fertilizer" = rgb(87, 62, 92, maxColorValue = 225),
# "+compost" = rgb(167, 156, 109, maxColorValue = 225)
# )
# )
custom_colors <- scale_fill_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
# Set a global theme
theme_set(theme_bw())
site.df <-
cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
mutate(
grazing_hist = ifelse(block == 1 | block == 2, "high", ifelse(block == 3 | block == 4, "low", NA))
) %>% # separate by block (grazing intensities)
group_by(nut_trt, ppt_trt, yr, code4, grazing_hist) %>% # make each code unique -- one value for species abundance for each trt
summarise(
# pct_cover,
abundance = mean(pct_cover)
# sd = sd(pct_cover)
)
`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.
matrix0 <- spread(
site.df,
code4,
abundance,
fill = 0)
# make rownames ----
rownames(matrix0) <- paste(matrix0$nut_trt, matrix0$ppt_trt, matrix0$yr, matrix0$grazing_hist, sep = "_")
Warning: Setting row names on a tibble is deprecated.
# delete columns ----
matrix <- matrix0[, c(5:79)]
# relativize data with Bray-Curtis dissimilarity
matrix.bray <- decostand(matrix, "total")
rownames(matrix.bray) <- rownames(matrix0)
Rules of thumb:
< 0.05 is excellent
0.05-0.1 is good
0.1-0.2 is fair
0.2-0.3 is cause for concern…
> 0.3 is poor, random
Stress: fair
0.1966553
Grouping by year and grazing history had most significant clustering shown.
nmds.comp <- metaMDS(matrix.bray, k=2, trymax = 25)
Run 0 stress 0.2020736
Run 1 stress 0.2023246
... Procrustes: rmse 0.07023633 max resid 0.190107
Run 2 stress 0.2017185
... New best solution
... Procrustes: rmse 0.02245434 max resid 0.1414201
Run 3 stress 0.2088616
Run 4 stress 0.1970656
... New best solution
... Procrustes: rmse 0.07244987 max resid 0.2329502
Run 5 stress 0.2176689
Run 6 stress 0.1970626
... New best solution
... Procrustes: rmse 0.006282932 max resid 0.03757378
Run 7 stress 0.2020248
Run 8 stress 0.2200448
Run 9 stress 0.2033604
Run 10 stress 0.1973254
... Procrustes: rmse 0.02289426 max resid 0.1158253
Run 11 stress 0.2042184
Run 12 stress 0.2020728
Run 13 stress 0.2017572
Run 14 stress 0.203447
Run 15 stress 0.2183276
Run 16 stress 0.1990614
Run 17 stress 0.2000129
Run 18 stress 0.2089696
Run 19 stress 0.2082397
Run 20 stress 0.2028091
Run 21 stress 0.2000005
Run 22 stress 0.1999666
Run 23 stress 0.2019172
Run 24 stress 0.2147543
Run 25 stress 0.2017483
*** Best solution was not repeated -- monoMDS stopping criteria:
25: stress ratio > sratmax
nmds.comp # pretty not great stress
Call:
metaMDS(comm = matrix.bray, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray
Distance: bray
Dimensions: 2
Stress: 0.1970626
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 6 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray’
stressplot(nmds.comp)
plot(nmds.comp, type="t")
nmds1 <- as.data.frame(scores(nmds.comp, choices=c(1), display=c("sites"))) %>%
mutate(matrix0[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.comp, choices=c(2), display=c("sites")))
nmds_dat <- cbind(nmds1, nmds2) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year")
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Grazing History")
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(grazing_hist), fill = as.factor(yr))) +
geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "Grouped by Grazing History & Year") +
scale_color_manual(
values = c(
"low" = "cyan",
"high" = "brown",
"2019" = "salmon",
"2020" = "limegreen",
"2021" = "cornflowerblue")
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(nut_trt)),
geom = "polygon",
fill = NA, # Makes the ellipse transparent inside
size = 1) +
labs(title = "Grouped by Nutrient Treatment & Year") +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
Stress: fair
0.1869068
# relativize data with Bray-Curtis dissimilarity
matrix.low <- matrix0 %>%
filter(grazing_hist == "low")
row.low <- paste(matrix.low$nut_trt, matrix.low$ppt_trt, matrix.low$yr, matrix.low$grazing_hist, sep = "_")
matrix.low <- matrix.low[, c(5:79)]
rownames(matrix.low) <- row.low
Warning: Setting row names on a tibble is deprecated.
matrix.bray.low <- decostand(matrix.low, "total")
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, color = as.factor(yr), shape = nut_trt)) +
geom_point(size = 3, stroke = 1) +
labs(title = "2019-2021 (LOW) - Grouping by Year & Nutrient Treatment")
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment")
# add year ellipses ----
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
custom_colors +
labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2019 <- matrix0 %>%
filter(yr == 2019)
row.2019 <- paste(matrix.2019$nut_trt, matrix.2019$ppt_trt, matrix.2019$yr, matrix.2019$grazing_hist, sep = "_")
matrix.2019 <- matrix.2019[, c(5:79)]
rownames(matrix.2019) <- row.2019
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2019 <- decostand(matrix.2019, "total")
Stress: Excellent
0.0560178
nmds.2019 <- metaMDS(matrix.bray.2019, k=2, trymax = 25)
Run 0 stress 0.0560178
Run 1 stress 0.05957995
Run 2 stress 0.05960972
Run 3 stress 0.05851815
Run 4 stress 0.05851812
Run 5 stress 0.05957995
Run 6 stress 0.06407773
Run 7 stress 0.05851805
Run 8 stress 0.05982901
Run 9 stress 0.06074817
Run 10 stress 0.05984339
Run 11 stress 0.05982901
Run 12 stress 0.07612032
Run 13 stress 0.06031567
Run 14 stress 0.05652852
Run 15 stress 0.06130516
Run 16 stress 0.05905532
Run 17 stress 0.0560178
... Procrustes: rmse 1.328847e-06 max resid 3.511506e-06
... Similar to previous best
Run 18 stress 0.07356079
Run 19 stress 0.05652852
Run 20 stress 0.05957995
*** Best solution repeated 1 times
nmds.2019 # pretty not great stress
Call:
metaMDS(comm = matrix.bray.2019, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2019
Distance: bray
Dimensions: 2
Stress: 0.0560178
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2019’
stressplot(nmds.2019)
plot(nmds.2019, type="t")
nmds1.2019 <- as.data.frame(scores(nmds.2019, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 4)
)
nmds2.2019 <- as.data.frame(scores(nmds.2019, choices=c(2), display=c("sites")))
nmds_2019_dat <- cbind(nmds1.2019, nmds2.2019) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
# read in trait data
# subset trait data to match comp data
# set matching row names
# relativize data with Bray-Curtis dissimilarity
matrix.2020 <- matrix0 %>%
filter(yr == 2020)
row.2020 <- paste(matrix.2020$nut_trt, matrix.2020$ppt_trt, matrix.2020$yr, matrix.2020$grazing_hist, sep = "_")
matrix.2020 <- matrix.2020[, c(5:79)]
rownames(matrix.2020) <- row.2020
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2020 <- decostand(matrix.2020, "total")
Stress: fair
0.1522495
nmds.2020 <- metaMDS(matrix.bray.2020, k=2, trymax = 25)
Run 0 stress 0.1524241
Run 1 stress 0.1729906
Run 2 stress 0.1595535
Run 3 stress 0.1584809
Run 4 stress 0.1524241
... New best solution
... Procrustes: rmse 1.001013e-05 max resid 2.563444e-05
... Similar to previous best
Run 5 stress 0.1522495
... New best solution
... Procrustes: rmse 0.02107748 max resid 0.06346541
Run 6 stress 0.1524241
... Procrustes: rmse 0.02107671 max resid 0.06351798
Run 7 stress 0.1522495
... New best solution
... Procrustes: rmse 6.786807e-06 max resid 1.943279e-05
... Similar to previous best
Run 8 stress 0.1821737
Run 9 stress 0.1594876
Run 10 stress 0.1872714
Run 11 stress 0.1705568
Run 12 stress 0.1752228
Run 13 stress 0.1522045
... New best solution
... Procrustes: rmse 0.01182953 max resid 0.03730308
Run 14 stress 0.1747777
Run 15 stress 0.1592892
Run 16 stress 0.1595535
Run 17 stress 0.1625903
Run 18 stress 0.1594876
Run 19 stress 0.165848
Run 20 stress 0.1594877
Run 21 stress 0.1524241
... Procrustes: rmse 0.0248137 max resid 0.06599096
Run 22 stress 0.1584809
Run 23 stress 0.1522495
... Procrustes: rmse 0.01183773 max resid 0.03738525
Run 24 stress 0.1744238
Run 25 stress 0.1744238
*** Best solution was not repeated -- monoMDS stopping criteria:
16: stress ratio > sratmax
9: scale factor of the gradient < sfgrmin
nmds.2020
Call:
metaMDS(comm = matrix.bray.2020, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2020
Distance: bray
Dimensions: 2
Stress: 0.1522045
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 13 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2020’
stressplot(nmds.2020)
plot(nmds.2020, type="t")
nmds1.2020 <- as.data.frame(scores(nmds.2020, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 4)
)
nmds2.2020 <- as.data.frame(scores(nmds.2020, choices=c(2), display=c("sites")))
nmds_2020_dat <- cbind(nmds1.2020, nmds2.2020) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 4, stroke = 1) +
labs(title = "2020 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 - Grouping by Nutrient Treatment")
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021 <- matrix0 %>%
filter(yr == 2021)
row.2021 <- paste(matrix.2021$nut_trt, matrix.2021$ppt_trt, matrix.2021$yr, matrix.2021$grazing_hist, sep = "_")
matrix.2021 <- matrix.2021[, c(5:79)]
rownames(matrix.2021) <- row.2021
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2021 <- decostand(matrix.2021, "total")
Stress: fair
0.1230831
nmds.2021 <- metaMDS(matrix.bray.2021, k=2, trymax = 25)
Run 0 stress 0.1232098
Run 1 stress 0.132972
Run 2 stress 0.132972
Run 3 stress 0.1506988
Run 4 stress 0.1515614
Run 5 stress 0.1232098
... New best solution
... Procrustes: rmse 1.844754e-06 max resid 4.667063e-06
... Similar to previous best
Run 6 stress 0.1230832
... New best solution
... Procrustes: rmse 0.04464562 max resid 0.1321501
Run 7 stress 0.1232098
... Procrustes: rmse 0.04464457 max resid 0.1315388
Run 8 stress 0.1421227
Run 9 stress 0.1339879
Run 10 stress 0.1515617
Run 11 stress 0.1421227
Run 12 stress 0.1506988
Run 13 stress 0.132972
Run 14 stress 0.1555173
Run 15 stress 0.1436877
Run 16 stress 0.1761755
Run 17 stress 0.1339879
Run 18 stress 0.1232098
... Procrustes: rmse 0.04464533 max resid 0.1315293
Run 19 stress 0.1507047
Run 20 stress 0.1525517
Run 21 stress 0.1339879
Run 22 stress 0.1230832
... Procrustes: rmse 2.173671e-05 max resid 6.138502e-05
... Similar to previous best
*** Best solution repeated 1 times
nmds.2021
Call:
metaMDS(comm = matrix.bray.2021, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2021
Distance: bray
Dimensions: 2
Stress: 0.1230832
Stress type 1, weak ties
Best solution was repeated 1 time in 22 tries
The best solution was from try 6 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2021’
stressplot(nmds.2021)
plot(nmds.2021, type="t")
nmds1.2021 <- as.data.frame(scores(nmds.2021, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 4)
)
nmds2.2021 <- as.data.frame(scores(nmds.2021, choices=c(2), display=c("sites")))
nmds_2021_dat <- cbind(nmds1.2021, nmds2.2021) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 - Grouping by Precipitation Treatment")
Greenhouse traits
Field traits
# read in trait data
# subset trait data to match comp data
# set matching row names